%--------------------------------------------------------------------------
% SIMULATION
%--------------------------------------------------------------------------

rng('default')
Nsim = 100; % number of firms


% simulate initial firm distribution without disasters,
% mimicking the 1947--2003 pre-estimation sample,
% based on which to initialize the actual simulation
burn_in         = 1000; % number of initial periods to discart
Tsim            = 10000;
[~,states_init] = sim_markov_chain_fast(P_nodis,Tsim+burn_in,1,1,1:2);
epsC_init       = randn(Tsim+burn_in,1);
epsCE_init      = [epsC_init,randn(size(epsC_init))]*chol([1,rho_CE;rho_CE,1]);
epsE_init       = epsCE_init(:,2);
g_sim           = mu_X(states_init) + sig_XS(states_init).*epsE_init + sig_Xid(states_init).*randn(Tsim+burn_in,Nsim); % [Tsim,Nsim]
k_sim_init      = NaN(size(g_sim));
k_sim_init(1,:) = k_bar(states_init(1));
k_QML_init      = NaN(size(g_sim));
k_QML_init(1,:) = k_bar(states_init(1));
last_issuance_init = ones(size(g_sim)); % times initial state of the Markov chain
for t = 1 : Tsim+burn_in-1 % maturity
    
    k_sim_init(t+1,:) = k_sim_init(t,:).*exp(g_sim(t,:));
    k_QML_init(t+1,:) = k_QML_init(t,:);
    
    default = k_sim_init(t,:) <= k_def(states_init(t));
    k_sim_init(t+1,default) = k_bar(states_init(t)).*exp(g_sim(t,default));
    k_QML_init(t+1,default) = k_bar(states_init(t));
    
    issuance = k_sim_init(t,:) >= k_iss(states_init(t));
    k_sim_init(t+1,issuance) = k_bar(states_init(t)).*exp(g_sim(t,issuance));
    k_QML_init(t+1,issuance) = k_bar(states_init(t));
    
    last_issuance_init(t:end,issuance)  = ones(Tsim+burn_in-t+1,sum(issuance)).*states_init(t);
    last_issuance_init(t:end,default)  = ones(Tsim+burn_in-t+1,sum(default)).*states_init(t);
    
end
k_sim_init(1:burn_in,:) = [];
k_QML_init(1:burn_in,:) = [];
last_issuance_init(1:burn_in,:) = [];
states_init(1:burn_in) = [];

% randomize cross-sections
ix = randperm(Tsim)';
k_sim_init = k_sim_init(ix,:);
k_QML_init = k_QML_init(ix,:);
last_issuance_init = last_issuance_init(ix,:);
states_init = states_init(ix);

% find periods corresponding to states 1 & 2
i1_init = find(states_init==1);
i2_init = find(states_init==2);


% simulate estimation sample
load aggregate_path_samples_long
delC = delC(:,1:5000);
epsC = epsC(:,1:5000);
states = states(:,1:5000);

Tsim   = size(states,1); % number of periods
Ssim   = size(states,2); % number of panels
epsCE  = [epsC(:),randn(size(epsC(:)))]*chol([1,rho_CE;rho_CE,1]); % account for correlation between consumption and systematic earnings shocks
epsE   = reshape(epsCE(:,2),size(epsC));

% simulate kappa panel
g_sim                 = mu_X(states) + sig_XS(states).*epsE + sig_Xid(states).*randn(Tsim,Ssim,Nsim); % [Tsim,Ssim,Nsim]
k_def_sim             = k_def(states); % [Tsim,Ssim]
k_iss_sim             = k_iss(states); % [Tsim,Ssim]
k_bar_sim             = k_bar(states); % [Tsim,Ssim]
i1                    = find(states(1,:)==1); % find periods corresponding to states 1 & 2
i2                    = find(states(1,:)==2);
k_sim                 = NaN(size(g_sim));
k_QML                 = NaN(size(g_sim));
last_issuance         = NaN(size(g_sim));
k_sim(1,i1,:)         = k_sim_init(i1_init(1:length(i1)),:);
k_sim(1,i2,:)         = k_sim_init(i2_init(1:length(i2)),:);
k_QML(1,i1,:)         = k_QML_init(i1_init(1:length(i1)),:);
k_QML(1,i2,:)         = k_QML_init(i2_init(1:length(i2)),:);
last_issuance(:,i1,:) = permute(repmat(last_issuance_init(i1_init(1:length(i1)),:),[1,1,Tsim]),[3,1,2]);
last_issuance(:,i2,:) = permute(repmat(last_issuance_init(i2_init(1:length(i2)),:),[1,1,Tsim]),[3,1,2]);
for t=1:Tsim-1
    
    k_sim(t+1,:,:)       = k_sim(t,:,:).*exp(g_sim(t,:,:));
    k_QML(t+1,:,:)       = k_QML(t,:,:);
    
    k_sim_temp           = k_sim(t+1,:,:); % temp: create 2D objects
    k_QML_temp           = k_QML(t+1,:,:);
    g_sim_temp           = g_sim(t,:,:);
    k_bar_temp           = k_bar(states(t,:))*ones(1,Nsim);
    states_temp          = states(t,:)'*ones(1,Nsim);
    
    default              = squeeze(k_sim(t,:,:)) <= k_def(states(t,:));
    k_sim_temp(default)  = k_bar_temp(default).*exp(g_sim_temp(default));
    k_QML_temp(default)  = k_bar_temp(default);
    issuance             = squeeze(k_sim(t,:,:)) >= k_iss(states(t,:));
    k_sim_temp(issuance) = k_bar_temp(issuance).*exp(g_sim_temp(issuance));
    k_QML_temp(issuance) = k_bar_temp(issuance);
    k_sim(t+1,:,:)       = k_sim_temp;
    k_QML(t+1,:,:)       = k_QML_temp;
    
    last_issuance_temp                = reshape(last_issuance(t:end,:,:),[Tsim-t+1,Ssim*Nsim]);
    last_issuance_temp(:,issuance(:)) = ones(Tsim-(t-1),1) * states_temp(issuance(:))';
    last_issuance(t:end,:,:)          = reshape(last_issuance_temp,[Tsim-t+1,Ssim,Nsim]);
    
end


default = k_sim <= repmat(k_def_sim,[1,1,Nsim]);


%--------------------------------------------------------------------------
% Table 10
%--------------------------------------------------------------------------


Fit     = griddedInterpolant({1:nS,k},lamDx);
DX_sim  = Fit(last_issuance,k_QML);

Fit     = griddedInterpolant({1:nS,k},lamSx);
SX_sim  = Fit(repmat(states,[1,1,Nsim]),k_sim);
SX_sim(default) = NaN;

Fit     = griddedInterpolant({1:nS,k},lamS);
S_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim);
S_sim(default) = NaN;

QMLEV_sim = DX_sim./( DX_sim + SX_sim .* k_sim./k_QML);
QMLEV_sim(default) = NaN;

Fit         = griddedInterpolant({1:nS,k},EP_R);
E_R_sim     = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},VP_R');
V_R_sim    = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k},IV_ATM);
IV_sim      = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit        = griddedInterpolant({1:nS,k},IV_Skew);
Skew_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim);

Fit         = griddedInterpolant({1:nS,k,1:nS},CDS(:,:,:,60));
CDS60_sim   = Fit(repmat(states,[1,1,Nsim]),k_sim,last_issuance);

ER_sim      = E_R_sim - reshape(Rf(states,1),size(states));

Rcum_sim    = exp(log(S_sim(2:end,:,:))  - log(SX_sim(1:end-1,:,:)) + g_sim(1:end-1,:,:)) - 1;

R_MKT_sim   = trimmean(Rcum_sim,2,3); % expected returns are measured contemporaneously to mktcaps

data = [ QMLEV_sim(:) ER_sim(:) CDS60_sim(:) IV_sim(:) Skew_sim(:) reshape(repmat([NaN(1,Ssim); R_MKT_sim],[1,1,Nsim]),[Tsim*Ssim*Nsim,1]) V_R_sim(:)];
del = any(isnan(data),2);
data(del,:) = [];

mean_vec    = mean(data(:,1:5));
var_vec     = var(data(:,1:6));
var_vec(2)  = mean(data(:,7)) + var(data(:,2));
table_10    = [mean_vec, sqrt(var_vec)]';
G_m         = [mean_vec, var_vec];

